!     Copyright (C) 2010 Benjamin Piaud
!
!     LIMBES is free software; you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation; either version 3, or (at your option)
!     any later version.

!     LIMBES is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.

!     You should have received a copy of the GNU General Public License
!     along with LIMBES; if not, see <http://www.gnu.org/licenses/>


!     subroutine LIMBES_write_macro write the macroscopic fields rho,ux,uy and T
!     in vtk legacy file format.
!     LIMBES_rho(1:LIMBES_nmax,1:LIMBES_nmax) -> rho (real*8)
!     LIMBES_ux(1:LIMBES_nmax,1:LIMBES_nmax) -> ux (real*8)
!     LIMBES_uy(1:LIMBES_nmax,1:LIMBES_nmax) -> uy (real*8)
!     LIMBES_T(1:LIMBES_nmax,1:LIMBES_nmax) -> T (real*8)
!     needed arguments
!     integer*8 time_index (index time of writen file)
!     integer*8 itime (simulation time)

subroutine LIMBES_write_macro(itime,time_index)
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  real*8 :: dx,Jx,Jy
  integer*8 :: itime,time_index,nz,xi,yi,ni
  integer*8 :: p1,p2,p3,p4,ncelltemp
  integer*8 :: nxmin,nxmax,nymin,nymax
  integer*8 :: num_zero
  character*50 :: file_namevtk
  character*6 :: czero,ctime

  nz=1
  dx=SIMUL_dx

  !     determine output region writing
  if (VTK_xmin.GT.0.OR.VTK_xmax.GT.0.OR.VTK_ymin.GT.0.OR.VTK_ymax.GT.0) then
     nxmin=VTK_xmin
     nxmax=VTK_xmax
     nymin=VTK_ymin
     nymax=VTK_ymax
     ncelltemp=0
     do yi=nymin,nymax
        do xi=nxmin,nxmax
           ni=LIMBES_cell_ni(xi,yi)
           if (ni.GT.0)  ncelltemp= ncelltemp+1
        enddo
     enddo
  else
     nxmin=1
     nxmax=SIMUL_nx
     nymin=1
     nymax=SIMUL_ny
     ncelltemp=SIMUL_cell_number
  endif

  czero='000000'
  write(ctime,'(I6)') time_index
  num_zero=intlen(SIMUL_save_number) - intlen(time_index)

  !     vtk format
  file_namevtk='macro'//czero(strlen(czero) -(num_zero-1):strlen(czero))//ctime(strlen(ctime)-(intlen(time_index)-1):strlen(ctime))//'.vtk'


  open(10,file=SIMUL_path(1:strlen(SIMUL_path))//file_namevtk(1:strlen(file_namevtk)))
  write(10,'(''# vtk DataFile Version 1.0'')')
  write(10,'(''LIMBES output'')')
  write(10,'(''ASCII'')')
  write(10,'(''           '')')
  write(10,'(''DATASET UNSTRUCTURED_GRID'')')
  write(10,'(''POINTS '',i10,'' float'')')(SIMUL_nx+1)*(SIMUL_ny+1)
  do yi=0,SIMUL_ny
     do xi=0,SIMUL_nx
        write(10,*) float(xi)*dx,float(yi)*dx,float(0)*dx
     enddo
  enddo
  write(10,'(''CELLS '',2i10)') ncelltemp,5*ncelltemp
  do yi=nymin,nymax
     do xi=nxmin,nxmax
        p1=(SIMUL_nx+1)*(yi-1) + (xi-1)
        p2=p1+1
        p3=p2 + SIMUL_nx + 1
        p4=p3 - 1
        if (LIMBES_cell_ni(xi,yi).GT.0) write(10,*) 4,p1,p2,p3,p4
     enddo
  enddo
  write(10,'(''CELL_TYPES '',i7)') ncelltemp
  do yi=nymin,nymax
     do xi=nxmin,nxmax
        ni=LIMBES_cell_ni(xi,yi)
        if (ni.GT.0)   write(10,'(i7)') 9
     enddo
  enddo
  write(10,'(''CELL_data '',i7)') ncelltemp
  write(10,'(''SCALARS rho float'')')
  write(10,'(''LOOKUP_TABLE default'')')
  do yi=nymin,nymax
     do xi=nxmin,nxmax
        ni=LIMBES_cell_ni(xi,yi)
        if (ni.GT.0) write(10,*) LIMBES_rho(ni)
     enddo
  enddo
  write(10,'(''VECTORS velocity float'')')
  do yi=nymin,nymax
     do xi=nxmin,nxmax
        ni=LIMBES_cell_ni(xi,yi)
        if (ni.GT.0) write(10,*) LIMBES_ux(ni),LIMBES_uy(ni),zero
     enddo
  enddo
  write(10,'(''SCALARS temperature float'')')
  write(10,'(''LOOKUP_TABLE default'')')
  do yi=nymin,nymax
     do xi=nxmin,nxmax
        ni=LIMBES_cell_ni(xi,yi)
        if (ni.GT.0) write(10,*) LIMBES_RT(ni)/FLUID_R
     enddo
  enddo
  write(10,'(''VECTORS heat_flux float'')')
  do yi=nymin,nymax
     do xi=nxmin,nxmax
        ni=LIMBES_cell_ni(xi,yi)
        if (ni.GT.0) then
           call heat_flux(ni,Jx,Jy)
           write(10,*) Jx,Jy,zero
        endif
     enddo
  enddo
  close(10)



end subroutine LIMBES_write_macro

!-------------------------------------------------------
!-------------------------------------------------------
!-------------------------------------------------------
!     subroutine LIMBES_save_f write the distribution function f in the file
!     LIMBES_f(1:LIMBES_nmax,1:LIMBES_nmax,1:LIMBES_nvmax,1:LIMBES_nvmax)
!     -> rho (real*8)

subroutine LIMBES_save_f()
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  integer*8 :: nv,vi,ni


  nv=LIMBES_nv
  
  open(10,file=SIMUL_path(1:strlen(SIMUL_path))//'f_save.dat',form='unformatted')       
  do ni=1,SIMUL_cell_number
     write(10) (LIMBES_f(ni,vi),vi=1,nv)
  enddo
  close(10)

end subroutine LIMBES_save_f
!-------------------------------------------------------
!-------------------------------------------------------
!-------------------------------------------------------
subroutine heat_flux(ni,Jx,Jy)
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  real*8 :: Jx,Jy,Pxy,vx,vy,ux,uy,w_f
  integer*8 :: vi,nv,ni

  nv=LIMBES_nv
  ux=LIMBES_ux(ni)
  uy=LIMBES_uy(ni)
  Jx=zero
  Jy=zero

  do vi=1,nv
     vx=LIMBES_vx(vi)
     vy=LIMBES_vy(vi)
     w_f=LIMBES_w(vi)*LIMBES_f(ni,vi)
     Jx=Jx + (vx-ux)*((vx-ux)**two + (vy-uy)**two)*w_f
     Jy=Jy + (vy-uy)*((vx-ux)**two + (vy-uy)**two)*w_f
     !     Jy=Jy + vy*(vx**two + vy**two)*w_f
     Pxy=Pxy +  (vx-ux)*(vy-uy)*w_f         
  enddo

  Jx=Jx/two
  Jy=Jy/two

end subroutine heat_flux
